#!/usr/bin/env python
#encoding: utf-8

import numpy as np
from sys import argv

from math import factorial
from math import gamma
import math

import pylab as pl

def at_least_1_of_z_out_of_x_total_when_y_is_taken(x,y,z):
#
	# Critical regions manually handled (math domain problems):
	
	x = int(x)
	y = int(y)
	z = int(z)
	total = 0

	for i in range(z): # i is "places of z taken"
		if x-z-y+i<0: continue # Continue - combination not possible
		if y-i<0: continue # Continue - combination not possible
		total = total + ( factorial(z) / (factorial(z-i)*factorial(i)) ) * ( factorial(x-z) ) / ( factorial((x-z-y)+i)*factorial(y-i) )
#
#		total = total + factorial(z) / factorial(z - i) / factorial(i) * factorial(x - z) / factorial(x - z - y + i) / factorial(y - i);

	# FUNCTION FOUND VIA MAPLE SIMPLIFY!
#	return -(-GAMMA(x + 1) * GAMMA(y + 1 - z) + GAMMA(x - z + 1) * GAMMA(y + 1)) / GAMMA(y + 1 - z) / GAMMA(x + 1);

	result = total / ( float(factorial(x)) / float( factorial(x-y)*factorial(y) ) )
#
	return result;


def create_transition_matrix(G,g,q,g1,g2):
	if g1+g2!=1: raise ValueError
	size = G+1
	P = np.matrix(np.zeros((size,size)))
	for row in range(size):
		if row == G:
			P[row,row] = 1;
		else:
#			print row, row
			P[row,row] = (g2/float(q**(G-row))) + (g1/float(q**((G-row)*at_least_1_of_z_out_of_x_total_when_y_is_taken(G,row,g) ) ) )
#			print row, row, row
			P[row+1,row] = 1 - P[row,row]
	return P
	
def plotting(P, no_packets):
	transfers = range(1,no_packets+1)
	s_0 = np.matrix(np.zeros( (P.shape[0],1) ) )
	s_0[0,0] = 1
	probability_of_full_rank = []
	
	for transfer in transfers:
		pmf = (P**transfer)*s_0
		pmf_full = pmf[-1,0]
		probability_of_full_rank.append(pmf_full)
		
#	print probability_of_full_rank
	
	pl.plot(transfers,probability_of_full_rank)

	
	
		

if __name__=="__main__":
#	print at_least_z_out_of_x_total_when_y_is_taken(64,64,64)
	if len(argv) > 3:
		print at_least_1_of_z_out_of_x_total_when_y_is_taken(argv[1],argv[2],argv[3])

	p2 = create_transition_matrix(40,15,256,1,0)
	plotting(p2, 60)
#	p1 = create_transition_matrix(15,15,256,0.1,0.9)
#	plotting(p1,200)
	
	pl.grid('on')
	pl.show()	
	
	
	
	
